In Silico Identification of a BRCA1:miR-29:DNMT3 Axis Involved in the Control of Hormone Receptors in BRCA1-Associated Breast Cancers

Germline inactivating mutations in the BRCA1 gene lead to an increased lifetime risk of ovarian and breast cancer (BC). Most BRCA1-associated BC are triple-negative tumors (TNBC), aggressive forms of BC characterized by a lack of expression of estrogen and progesterone hormone receptors (HR) and HER2. How BRCA1 inactivation may favor the development of such a specific BC phenotype remains to be elucidated. To address this question, we focused on the role of miRNAs and their networks in mediating BRCA1 functions. miRNA, mRNA, and methylation data were retrieved from the BRCA cohort of the TCGA project. The cohort was divided into a discovery set (Hi-TCGA) and a validation set (GA-TCGA) based on the platform used for miRNA analyses. The METABRIC, GSE81002, and GSE59248 studies were used as additional validation data sets. BCs were differentiated into BRCA1-like and non-BRCA1-like based on an established signature of BRCA1 pathway inactivation. Differential expression of miRNAs, gene enrichment analysis, functional annotation, and methylation correlation analyses were performed. The miRNAs downregulated in BRCA1-associated BC were identified by comparing the miRNome of BRCA1-like with non-BRCA1-like tumors from the Hi-TCGA discovery cohort. miRNAs:gene-target anticorrelation analyses were then performed. The target genes of miRNAs downregulated in the Hi-TCGA series were enriched in the BRCA1-like tumors from the GA-TCGA and METABRIC validation data sets. Functional annotation of these genes revealed an over-representation of several biological processes ascribable to BRCA1 activity. The enrichment of genes related to DNA methylation was particularly intriguing, as this is an aspect of BRCA1 functions that has been poorly explored. We then focused on the miR-29:DNA methyltransferase network and showed that the miR-29 family, which was downregulated in BRCA1-like tumors, was associated with poor prognosis in these BCs and inversely correlated with the expression of the DNA methyltransferases DNMT3A and DNMT3B. This, in turn, correlated with the methylation extent of the promoter of HR genes. These results suggest that BRCA1 may control the expression of HR via a miR-29:DNMT3:HR axis and that disruption of this network may contribute to the receptor negative phenotype of tumors with dysfunctional BRCA1.


Introduction
Patients carrying inactivating BRCA1 germline mutations have an increased lifetime risk of developing breast and ovarian cancer [1]. Specifically, BRCA1 germline mutations predispose to the occurrence of triple-negative breast cancer (TNBC), an aggressive subtype of breast cancer (BC) that is negative for the expression of the hormone receptors (HR), ERα and PR, and HER2 [2,3].
In addition to germline mutations, which account for approximately 5-10% of all BCs, a subset of sporadic cancers, most of the TNBC subtype, have features of dysfunctional BRCA1 [4][5][6]. Approximately 10% exhibit aberrant methylation of the BRCA1 promoter, a Int. J. Mol. Sci. 2023, 24, 9916 2 of 18 feature closely associated with the downregulation of BRCA1 gene transcripts and loss of heterozygosity at the BRCA1 locus [7][8][9]. This fraction of tumors shows hallmarks of dysfunctional BRCA1 characteristics of hereditary BRCA1 tumors, such as the absence of HR, amplification of MYC, and an unchanged ERBB2 gene (encoding the HER2 receptor) [10,11]. Finally, patterns of genomic alterations typical of familial BRCA1 tumors are detected in an additional fraction of sporadic BCs devoid of evident BRCA1 gene or expression alterations [11,12]. Thus, in addition to hereditary BRCA1 tumors, a substantial proportion of sporadic BCs may be considered BRCA1-associated since they rely on dysfunctional BRCA1 as a driving force. To better identify these tumors, Chen and coworkers developed a specific signature that they used to define a group of BRCA1-like tumors [5].
The large nuclear BRCA1 protein acts as a guardian of chromosome integrity through various functions that ensure the assembly and activity of macromolecular complexes involved in DNA double-strand break repair [13] and mitotic and replication control [14]. Tumors with deleterious mutations in the BRCA1 gene and hence impaired DNA repair capacity are highly sensitive to inter-strand crosslinking agents (platinum or alkylating agents), topoisomerase II inhibitors (anthracyclines), and PARP inhibitors [15][16][17].
Increasing evidence indicates that BRCA1 also plays a central role in regulating the transcription and processing of RNA [18,19], including non-coding microRNAs (miR-NAs) [20]. miRNAs are important regulators of gene expression by inhibiting the translation of target mRNAs, mostly by binding to their 3 UTR region [21,22]. BRCA1 is involved in both the induction and repression of miRNA transcription and maturation through the binding to DROSHA [20,[23][24][25]. Indeed, several miRNAs have been identified as deregulated in BRCA1/2-associated BCs [26,27], but how the interplay between BRCA1 and miRNAs may promote the development of the BRCA1-like/TNBC subtype remains unclear.
To identify miRNA:mRNA circuits involved in the induction of a BRCA1-like phenotype, we interrogated publicly available transcriptional profiling data. This approach allowed us to unveil a particular interplay between the miR-29 family of miRNAs, DNA methyltransferases, and HR.

Downregulated miRNAs in BRCA1-like Tumors
To address the hypothesis of a role for BRCA1:miRNA networks in the development of the BRCA1-like/TNBC subtype, microRNA expression data were retrieved from the TCGA-BRCA project (https://www.cancer.gov/tcga) [28]. To capture all cases with a putative dysfunctional BRCA1, we used the classification into BRCA1-like and non-BRCA1-like cases proposed by Chen and coworkers [5]. We used the TCGA BC miRNA isoform data generated by Hiseq Illumina technology (Hi-TCGA data set) as the discovery data set. This consists of 554 patients, 161 BRCA1-like, and 393 non-BRCA1-like. Differential expression analysis showed that 57 miRNAs were downregulated and 105 were upregulated in BRCA1like compared with non-BRCA1-like samples (log2FC > |0.4|, padj < 0.00001; Figure 1A and Supplemental Table S1).
To uncover genes whose overexpression, as a consequence of the alleviation from miRNAs-mediated inhibition, may promote tumorigenic properties, we focused on downregulated miRNAs and queried the miRDIP platform (version 5.0.2.3, https://ophid. utoronto.ca/mirDIP/; accessed on 20 April 2022). This is a large database for predicting miRNA:gene target pairs based on 3 UTR pairing [29]. The 57 downregulated miRNAs were predicted to target 4236 genes with high confidence (very high score). We assumed that if a miRNA acts as a negative regulator of a given mRNA, its expression levels should be somehow anticorrelated. Therefore, we retrieved the RNAseq data for the same cases of the Hi-TCGA BC series and performed anticorrelation analyses with miRNA expression. Using a correlation threshold R < −0.35 as a cutoff, 2215 interactions were observed between 41 of the 57 miRNAs used as input and 504 genes (Supplemental Table S2). Some genes were targeted by nearly 30 miRNAs ( Figure 1B). After having identified these anticorrelated miRNA:gene pairs using the Hi-TCGA discovery cohort, we sought to verify whether the genes potentially targeted by the miRNAs downregulated in BRCA1-like cases were enriched in BRCA1-like BC. We then performed gene set enrichment analysis (GSEA) on the mRNA data obtained from two independent data sets: the first, consisting of 76 BRCA1-like and 149 non-BRCA1-like BCs, was generated by the TCGA consortium using a genome analyzer sequencing platform to produce miRNA data (GA-TCGA); the second data set was obtained from the METABRIC project and included 215 BRCA1-like and 1068 non-BRCA1-like BCs [30]. The GSEA revealed that the 504 genes potentially targeted by downregulated miRNAs were indeed significantly enriched in BRCA1-like samples of both validation series ( Figure 1C,D), supporting the hypothesis of a BRCA1-mediated control over these genes via miRNAs.
Functional annotation by overrepresentation analysis (ORA) of the genes targeted by the downregulated miRNAs revealed that many of the over-represented biological processes of gene ontology (GO) were known to be controlled by BRCA1, such as cell cycle and cell division and DNA metabolic processes. Among the genes belonging to these processes, we noticed the involvement of the DNA methyltransferases, including DNMT3A and DNMT3B ( Figure 1E, Supplemental Table S3; Supplemental Figure S1). The two DNA methyltransferases catalyze the methylation of genomic DNA to establish new DNA methylation patterns during embryogenesis [31] and in cancer [32,33].

miR-29s:DNMT3A-DNMT3B Network Is Involved in BRCA1-like Tumors
Correlation analysis revealed an inverse relationship between DNMT3A and DNMT3B and a limited number of miRNAs downregulated in BRCA1-like tumors. The miR-29 family was strongly represented (Figure 2A; Supplemental Table S2) and capable of affecting both DNMT3A and DNMT3B. Indeed, the anti-correlation between miR-29b-2-5p, miR-29c-3p, and miR-29c-5p transcripts and DNMT3A and DNMT3B levels was demonstrated not only in the Hi-TCGA discovery series but also in the GA-TCGA validation data set and in two other cohorts containing data for both miRNAs and coding mRNA, GSE59248 (39 patients, [34]) and GSE81002 (377 patients, [35]), ( Figure 2B). DNMT3A did not result in being anticorrelated with the three miRNAs in the METABRIC series; nevertheless, a high anti-correlation between miR-29s and DNMT3B was confirmed in this cohort ( Figure 2B). Therefore, we investigated the role of miR-29s in BC and performed a more detailed analysis of the miR-29s:DNMT3A-DNMT3B network. First, we confirmed the downregulation of miR-29b-2-5p, miR-29c-3p, and miR-29c-5p in BRCA1-like BCs in both GA-TCGA and METABRIC validation sets ( Figure 3A,B).
We then examined the prognostic role of miR-29 members in both TCGA and METABRIC cohorts. Due to the limited number of BRCA1-like samples with high miR-29s levels, we pooled Hi-TCGA and GA-TCGA data after having calculated miR-29s z-scores. Intriguingly, the BRCA1-like group with low levels of the three miR-29s demonstrated a trend toward the worst disease-free survival (DFS) probability in both the TCGA ( Figure 3C) and METABRIC ( Figure 3D) data sets. The differences between groups reached statistical significance in the case of miR-29c-3p, the most highly expressed of the three miRNAs. The clinical relevance of miR-29s may involve DNMT3A and DNMT3B, which were indeed upregulated in both the Hi-TCGA discovery group and the GA-TCGA validation cohort ( Figure 4A,B). Consistent with the correlation data, in the METABRIC series, only DNMT3B was significantly increased in BRCA1-like cancers compared with non-BRCA1-like samples ( Figure 4C). series; nevertheless, a high anti-correlation between miR-29s and DNMT3B was confirmed in this cohort ( Figure 2B). Therefore, we investigated the role of miR-29s in BC and performed a more detailed analysis of the miR-29s:DNMT3A-DNMT3B network. First, we confirmed the downregulation of miR-29b-2-5p, miR-29c-3p, and miR-29c-5p in BRCA1-like BCs in both GA-TCGA and METABRIC validation sets ( Figure 3A,B). Interestingly, DNMT3s expression is associated with typical features of BRCA1-like samples, such as negativity for HR. Indeed, DNMT3B transcripts were significantly and markedly higher in samples classified as negative for either estrogen receptor-α (ERα) or progesterone receptor (PR) in both the Hi-TCGA discovery cohort and the GA-TCGA and METABRIC validation series ( Figure 5A,F, lower panels). DNMT3A transcripts were significantly higher in hormone receptor-negative samples from both Hi-TCGA and GA-TCGA ( Figure 5A,B,D,E, upper panels), but not in the METABRIC series ( Figure 5C,F), consistent with the lack of correlation with miR-29 family members and the absence of upregulation in BRCA1-like tumors in this cohort. We further examined the GSE81002 and GSE59248 series: DNMT3A and DNMT3B levels were higher in ERα-negative than in positive cases in both series (Supplemental Figures S1 and S2A) and in PR-negative samples in GSE59248 (Supplemental Figure S1B; the GSE81002 project did not provide PR status).   Samples were further split into BRCA1-like (orange lines) and non-BRCA1-like (blue lines) samples. A log-rank test was used to determine the differences between curves. The numbers under the graphs represent patients at risk at a specific time point (months). Survival curves are truncated at 100 months. and METABRIC ( Figure 3D) data sets. The differences between groups reached statistical significance in the case of miR-29c-3p, the most highly expressed of the three miRNAs.
The clinical relevance of miR-29s may involve DNMT3A and DNMT3B, which were indeed upregulated in both the Hi-TCGA discovery group and the GA-TCGA validation cohort ( Figure 4A,B). Consistent with the correlation data, in the METABRIC series, only DNMT3B was significantly increased in BRCA1-like cancers compared with non-BRCA1like samples ( Figure 4C). Interestingly, DNMT3s expression is associated with typical features of BRCA1-like samples, such as negativity for HR. Indeed, DNMT3B transcripts were significantly and markedly higher in samples classified as negative for either estrogen receptor-α (ERα) or progesterone receptor (PR) in both the Hi-TCGA discovery cohort and the GA-TCGA and METABRIC validation series ( Figure 5A,F, lower panels). DNMT3A transcripts were significantly higher in hormone receptor-negative samples from both Hi-TCGA and GA-TCGA ( Figure 5A,B,D,E, upper panels), but not in the METABRIC series ( Figure 5C,F), consistent with the lack of correlation with miR-29 family members and the absence of upregulation in BRCA1-like tumors in this cohort. We further examined the GSE81002 and GSE59248 series: DNMT3A and DNMT3B levels were higher in ERα-negative than in Statistical differences were determined with the Welch two-sample t-test.

DNA Methyltransferases Affect ESR1 and PGR DNA Methylation in BRCA1-like Samples
The data reported so far suggest a potential form of control over DNMT3s by BRCA1 via miR-29 family members and the association of DNMT3s expression with the absence of HR.
Interestingly, genes encoding ERα (ESR1 gene) and PR (PGR gene), which are weakly expressed in BRCA1-associated tumors [36], have been reported to be downregulated by methylation of promoter CpGs [37][38][39][40]. Therefore, we sought to determine whether the upregulation of DNA methyltransferases was associated with promoter methylation of ESR1 and PGR genes. To investigate this hypothesis, we obtained from the TCGA-BRCA project the methylation data of the genomic regions where ESR1 and PGR were mapped ( Figure 6A,B). The ESR1 gene spans more than 400 Kbp, includes 8 coding exons and, through the use of alternative promoters and splicing, forms dozens of transcript variants (RefSeq: NCBI Reference Sequence Database), which in turn are translated into at least 4 isoforms [41]. The PGR gene encodes two major isoforms, isoform A (PRA) and B (PRB), which are derived from two alternative promoters and two different transcription start sites (TSS) [42].
The TCGA methylation data consisted of β-values, the estimates of methylation levels, for 63 probes along the ESR1 gene and 17 probes for the PGR gene.
To determine which methylated CpGs might be relevant to control ESR1 and PGR transcript levels, we used data from Hi-TCGA patients (433 samples) and categorized hormone-positive and -negative samples based on the pit of the bimodal distribution of ESR1 and PGR log2RSEM-normalized counts (Supplemental Figure S3). A significant fold change >1.5 of CpG β-values in HR negative over positive samples was used to select CpGs whose methylation inhibits ESR1 and PGR expression. A limited number of CpG probes (Supplemental Table S4) met these criteria for the ESR1 gene: seven CpG probes (CpGs_A; Figure 6A) located near the transcription start sites (TSSs) of RefSeq NM_001122740 and NM_001385570 (variants 4 and 10, respectively) located far upstream of the canonical NM_000125 variant 1; a CpG (CpG_B) in the region of promoter B [38] near the TSSs of NM_001122740 and NM_001385569 (variants 2 and 9, respectively) and a CpG (CpG_C) close to TSS of NM_001328100 (variant 7). The available data for the CpG island were very limited and included only 5 of the 105 CpGs reported in the CpG islands track of the UCSC genome browser (https://genome.ucsc.edu/; accessed on 1 December 2022). These CpGs were only weakly methylated and appeared to scarcely affect ESR1 expression in the BC series of TCGA. positive cases in both series (Supplemental Figures S1 and S2A) and in PR-negative samples in GSE59248 (Supplemental Figure S1B; the GSE81002 project did not provide PR status).  project the methylation data of the genomic regions where ESR1 and PGR were mappe ( Figure 6A,B). The ESR1 gene spans more than 400 Kbp, includes 8 coding exons an through the use of alternative promoters and splicing, forms dozens of transcript varian (RefSeq: NCBI Reference Sequence Database), which in turn are translated into at least isoforms [41]. The PGR gene encodes two major isoforms, isoform A (PRA) and B (PRB which are derived from two alternative promoters and two different transcription sta sites (TSS) [42].  In the case of the PGR gene, two CpGs tested by TCGA probes (Supplemental Table  S5), CpG_A and CpG_B, showed β-values FC > 1.5 and were selected for subsequent analyses ( Figure 6B). In addition, the methylation of CpG islands seems to control PGR expression, even though the number of probes was limited (3 of 32 for CpGs-island 32 and 4 of 170 for CpGs-island170).
Next, we wanted to test whether DNMT3A and DNMT3B might have a role in the methylation of the selected CpGs in the ESR1 and PGR regions. To this end, we examined the correlation between the average β-values of the 7 CpGs in CpGs_A, the single CpG_B, and CpG_C β-values of ESR1 and the expression levels of DNMT3A and DNMT3B. Only CpG_C significantly correlated with DNMT3A in both the Hi-TCGA ( Figure 7A) and GA-TCGA (90 samples; Figure 7B) series. The low levels of CpG methylation of the CpG island did not correlate with DNMT3s expression in both Hi-TCGA and GA-TCGA cohorts (Supplemental Figure S4).    Supplemental Tables S4 and S5. Interestingly, both CpG_A and CpG_B of the PGR gene correlated markedly and significantly with both DNMT3A and DNMT3B in the Hi-TCGA and GA-TCGA cohorts ( Figure 7C,D). The levels of DNMT3A and DNMT3B transcripts were not associated with the methylation of CpG islands, which were barely methylated in either gene ( Figure 7C,D).
Thus, DNMT3A and DNMT3B were associated with the methylation of a number of CpGs that control the expression of HR genes, particularly of PGR.
Finally, to come full circle, we showed that the methylation levels of CpGs_C in the ESR1 gene and CpG_A and CpG_B in the PGR gene were indeed higher in BRCA1-like compared with BRCA1-non-like tumors in both Hi-TCGA and GA-TCGA series ( Figure 8A,B). Interestingly, both CpG_A and CpG_B of the PGR gene correlated markedly and significantly with both DNMT3A and DNMT3B in the Hi-TCGA and GA-TCGA cohorts (Figure 7C,D). The levels of DNMT3A and DNMT3B transcripts were not associated with the methylation of CpG islands, which were barely methylated in either gene ( Figure 7C,D).
Thus, DNMT3A and DNMT3B were associated with the methylation of a number of CpGs that control the expression of HR genes, particularly of PGR.
Finally, to come full circle, we showed that the methylation levels of CpGs_C in the ESR1 gene and CpG_A and CpG_B in the PGR gene were indeed higher in BRCA1-like compared with BRCA1-non-like tumors in both Hi-TCGA and GA-TCGA series ( Figure  8A,B).

Discussion
Inheritance of a germline mutation of the BRCA1 gene predisposes with high penetrance to breast, ovarian, and other forms of epithelial carcinoma [43]. Interestingly, the vast majority of breast cancers that occur in BRCA1 mutation carriers have a TNBC phenotype that prevents patients from being treated with therapies targeting hormone receptor pathways [3]. Despite much progress in elucidating how dysfunctional BRCA1 causes the occurrence of such specific BCs, we still do not have a complete picture. It is known that BRCA1 can both inhibit ERα signaling and induce ESR1 gene expression by binding coactivators such as BRD7 and OCT1 or resolving R-loops in the enhancer region upstream of the ESR1 gene, allowing ESR1 transcription [44][45][46][47][48]. In addition, it can regulate the PR pathway indirectly via ESR1, which induces transcription of the PGR gene, and also interact directly with PRA and PRB [42,49,50].

Discussion
Inheritance of a germline mutation of the BRCA1 gene predisposes with high penetrance to breast, ovarian, and other forms of epithelial carcinoma [43]. Interestingly, the vast majority of breast cancers that occur in BRCA1 mutation carriers have a TNBC phenotype that prevents patients from being treated with therapies targeting hormone receptor pathways [3]. Despite much progress in elucidating how dysfunctional BRCA1 causes the occurrence of such specific BCs, we still do not have a complete picture. It is known that BRCA1 can both inhibit ERα signaling and induce ESR1 gene expression by binding coactivators such as BRD7 and OCT1 or resolving R-loops in the enhancer region upstream of the ESR1 gene, allowing ESR1 transcription [44][45][46][47][48]. In addition, it can regulate the PR pathway indirectly via ESR1, which induces transcription of the PGR gene, and also interact directly with PRA and PRB [42,49,50].
Here, we propose that the downregulation of specific miRNAs, such as miR-29 family members, contributes to the lack of HR expression in BRCA1-like tumors.
There is increasing evidence that BRCA1 controls the expression of miRNAs through various mechanisms: transcriptional repression, activation, or by controlling miRNA biogenesis via the DROSHA microprocessor complex [20,23,25]. Either way, in the present analyses, we observed several miRNAs deregulated in BRCA1-like tumors. The targets of these miRNAs were involved in pathways known to be influenced by BRCA1 mutations, namely cell cycle, cell division, and DNA replication. Of note, the DNA methyltransferases DNMT3A and DNMT3B emerged among the genes involved in these pathways. They anticorrelated with miR-29 family members and were upregulated in BRCA1-like cancers. In vitro inhibition of these methyltransferases by miR-29 family members has already been demonstrated for several pathologies [51,52].
Although we could not demonstrate in these in silico analyses whether and how BRCA1 directly regulates miR-29 expression, a recent paper reported that BRCA1 binds cis-regulatory elements in the promoter region of miR-29b-1 (in a cluster with miR-29a on chromosome 7) and induces its expression [53]. We did not verify whether this form of control also applies to the other miR-29s, but it is worth noting that the genes encoding the most downregulated miR-29c -3p and -5p, and miR-29b-2-5p are clustered together on chromosome 1, suggesting a common regulation for these miRNAs. In this scenario, BRCA1 dysfunction could be the cause of low levels of these miR-29 family members. Interestingly, several lines of evidence suggest that miR-29s are downregulated in cancer and that their reduction is associated with a poor prognosis [54][55][56]. Here, we showed that the prognostic role of miR-29s was mainly associated with BRCA1-like samples. Indeed, low levels of these miRNAs, especially miR-29c-3p, in the BRCA1-like samples showed the worst prognosis and were associated with the short-term relapse that is typical of TNBC.
DNMT3A and DNMT3B are de novo methyltransferases that transfer a methyl group to the C-5 position of the cytosine residue to establish DNA methylation [61,62]. They may control gene expression by methylating their promoters, a mechanism that has long been considered critical for gene silencing in cancer development.
Even the MIR29C and MIR29B2 regions were methylated but the extent of methylation did not correlate with either miRNA expression or with DNMTs. Curiously, methylation of CpGs 20kbp upstream of miRNAs, which Poli E et al. reported to be inversely correlated with miR-29c and miR-29b-2 expression [63], actually partially correlated with the poorly expressed strands (mir-29c-5p and miR-29b-2-5p). No correlation was observed between the methylation of CpGs in the miR-29a/miR-29b-1 locus and miRNA expression. Thus, overall, our analyses suggest that the anti-correlation between DNMT3A/3B and miR-29s is not a consequence of DNMT3A/B controlling miR-29s expression but rather the opposite. Nevertheless, we cannot rule out a negative feedback loop in which DNMTs might partially control miR-29c expression by methylating its gene.
Remarkably, the data reported here showed a correlation between DNMT3A and DNMT3B expression and negativity for ERα and PR and a direct correlation with promoter methylation levels of HR, particularly the PGR gene, suggesting a BRCA1-mediated control of these genes via the miR-29s/DNMT3A-DNMT3B axis.
In the ESR1 gene, the CpG most affected by the miR-29s/DNMT3A-DNMT3B circuit is CpG_C, which maps close to the TSS of transcripts that produce the 46 and 36 kDa ERα truncated isoforms (ERα46 and ERα36, respectively). Both isoforms have been reported to be expressed in BC, even in TNBC, and to exert opposite functions in controlling cancer progression: the best studied ERα36 is capable of activating multiple signaling pathways critical for cancer aggressiveness and metastatic potential, whereas ERα46 associates with lower grade and smaller tumors [64][65][66][67]. Further studies are needed to characterize the isoforms affected by CpG DNA methylation and to discover the role of promoter methylation/expression of these isoforms in BC progression.
In the case of PGR, DNMT3A and DNMT3B are strongly associated with methylation of the two CpGs that appear to inhibit PGR expression, even if one of them maps on an exon far from the TSS of the PGR gene. Several lines of evidence suggest that CpG methylation plays a role in controlling gene expression even when far from the promoter by altering chromatin structure and binding of transcription factors [68,69]. Importantly, the two CpGs in the PGR gene were more methylated in BRCA1-like tumors than in non-BRCA1-like tumors.
We recognize that this work has the limitations of an in silico study that showed correlations and, thus, needs further in vitro validation to confirm the results. In addition, we used the BRCA1-like and non-BRCA1-like classification as a surrogate of a categorization that considered all types of BRCA1 inactivation (e.g., mutation, methylation, small or large deletion). The fact that miR-155, a miRNA known to be repressed by BRCA1, emerged among the miRNAs overexpressed in BRCA1-like supports the efficacy of this categorization. Thus, by analyzing a series of data sets, this study identifies a potential network in which the absence of a functional BRCA1 epigenetically suppresses HR. This network is thus implicated in the development of the typical BRCA1-associated tumors that barely express ERα and PR. Furthermore, this study extends the knowledge of DNA methylation control over ESR1 gene expression in the context of BRCA1-like cancers, as reported for a limited set of CpGs in the seminal work of Archey and coworkers [70]. Interestingly, promoter methylation and lack of HR expression lead to poor prognosis and resistance to a number of anti-hormonal therapies in BC [38,40,71]. Therefore, the network described here, in which dysfunctional BRCA1 leads to the downregulation of miR-29 accompanied by an increase in DNMT3A and DNMT3B and hormone gene methylation, paves the way for the development of novel strategies to overcome endocrine resistance in BCs.

Data Sets
In selecting data sets for this study, we first chose those that contained mRNA and miRNA data and the BRCA1 status of the same samples. Additional data sets, even without BRCA1 status information, were selected for specific validation (see the workflow in Supplemental Figure S5).
For miRNA isoforms, we used both the read counts from the Hi-TCGA data set and the normalized count in reads-per-million-miRNA-mapped from the Hi-TCGA and GA-TCGA data sets (TCGA miRNA pipeline at https://docs.gdc.cancer.gov/Data/Bioinformatics_ Pipelines/miRNA_Pipeline/). For survival analyses, we used both the Hi-TCGA and GA-TCGA data sets. We calculated the z-score of the log2 of the normalized count in reads-per-million-miRNA-mapped for each cohort: z-score = (miR-X sampleA − mean miR-X allSamples)/ds miR-X allSamples Subsequently, the z-scores of the two platforms were merged. Gene-level transcription estimates from the RNA-seq data were expressed as log2 RSEM-normalized counts. DNA methylation data report the methylation score for each CpG site expressed as β-values, which ranged from zero to one, corresponding to unmethylated and fully methylated DNA, respectively (https://docs.gdc.cancer.gov/Data/ Bioinformatics_Pipelines/). We analyzed the β-values of 63 probes in the ESR1 gene region (chr6: 152,011,103-152421432; GRCh37/hg19) and of 17 probes in the PGR region (chr11: 100,905,070-101001533; GRCh37/hg19).

Other Validation Data Sets
Two independent data cohorts with miRNA, mRNA, and clinicopathological data were identified in GEO data sets (https://www.ncbi.nlm.nih.gov/gds; accessed on 18 May 2022): the GSE81002 and GSE59248 series. GSE81002 included data from the Human miRNA Microarray v14 Rev.2 and the Agilent-028004 SurePrint G3 Human GE 8 × 60K Microarray for 377 patients [35]. GSE59248 contained data from the same platform as GSE81002 for mRNA quantification and data from Agilent-031181 Unre-stricted_Human_miRNA_V16.0_Microarray 030, 840, for 39 samples [34]. Since these data sets contained data for multiple probe sets, we collapsed probes into miRNA by selecting the probe with the highest average expression [73].

Differential Expression Analysis
We performed differential expression analysis using the Hi-TCGA data set, which includes 161 BRCA1-like and 393 non-BRCA1-like samples. Read counts of miRNA isoforms were used for the analysis, which was performed using the DESeq2 R package version 1.36.0 [74]. miRNAs with log2FC > |0.4|, BH-adjusted p values (padj) <0.00001, and base mean >5 normalized counts were considered differentially expressed.
Gene ontology (GO) enrichment analysis for the GO Biological Processes was performed by taking advantage of the Metascape tool (https://metascape.org; [77]). Pathway and process enrichment was performed.

Disease-Free Survival and Statistical Analysis
Spearman rank order correlation coefficients were computed between the normalized miRNA expression levels from the Hi-TCGA data set and the corresponding log2 RSEMnormalized counts of the putative targets. Statistical differences between groups were evaluated using the Welch t-test or a Mann-Whitney rank sum test, depending on the normal distribution of the variables considered. Tests were performed with R functions (R-3.6.2). Disease-free survival analyses were conducted considering updated clinical data (disease-free months and events) obtained from the cBioPortal (https://www.cbioportal. org/; accessed on 10 May 2023) for both TCGA and METABRIC cohorts. Patients were classified according to BRCA1 status and high or low levels of miR-29s, which were in turn determined by calculating the best cut-point based on the maximum log-rank value using the R package Survminer [78]. Survival was estimated using the Kaplan-Meier method, and differences between groups were evaluated with the log-rank test using the Survival and Survminer R packages [78,79].